C.................... CMINOR.FOR;16        26-APR-1994 14:53:26.59 
C....... This subroutine evaluates the minor ion and neutral densities
C....... from chemical equilibrium and is also used for printing the
C....... production and loss rates for the FLIP model. P. Richards
C....... Sep 1989
      SUBROUTINE CMINOR(I,J,JP,IPR,FD,ID,HE)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      DIMENSION PROD(6),LOSS(6),RTS(199),FEX(2),FD(9)
C----- commons to hold the EUV, e* and AURoral production rates
      INTEGER IVERT
      REAL EUVION,PEXCIT,PEPION,PAUION,PAUEXC,SUMION,SUMEXC
      COMMON/EUVPRD/EUVION(3,12,401),PEXCIT(3,12,401),PEPION(3,12,401)
      COMMON/BPRAUR/PAUION(3,12,401),PAUEXC(3,12,401)
      COMMON/BPRTOT/SUMION(3,12,401),SUMEXC(3,12,401)
      COMMON/MIA/ZNE(401),NOP(401),N2PLS(401),O2P(401),NPLUS(401)
      COMMON/MIB/OP2D(401),OP2P(401),N2D(401),N4S(401),NNO(401)
      COMMON/VODDN/VN2D(401),VN4S(401),VNNO(401),VO1D(401),IVERT
      COMMON/MIC/O1D(401),N2P(401)
      COMMON/FJS/N(4,401),TI(3,401),F(20)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/LPS/SZA(401),EPSN,DC,ISPC,I1,I2,IPV,JTIMAX,IPP,ISKP
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      COMMON/SAV/NSAVE(2,401),TISAV(3,401),FY(2,401),UN(401),EHT(3,401)
      COMMON/RK/VC(401)
      COMMON/NPRD/COLUM(3,401),OTHPR1(6,401),OTHPR2(6,401),OTHPR3(6,401)
      COMMON/MINOR/XIONN(9,401),XIONV(9,401),XMASS(9)
C....... common suggested by Kent Miller
      COMMON/CMCOM/PO2P
      DATA DNOP,DO2P,HEPLUS/3*0.0D0/
C
C....... for the production rates consult subroutine prdint
C
      DO 44 II=1,9
 44   FD(II)=0.0D0
      IF(Z(J).GT.900) RETURN
      TOP2PI=SUMION(1,9,J)
      TOP2DI=SUMION(1,8,J)
      TOTO2I=SUMION(2,7,J)
      DISNP=SUMION(3,4,J)+SUMION(3,5,J)+SUMION(3,6,J)
      PEDISNP=PEPION(3,4,J)+PEPION(3,5,J)+PEPION(3,6,J)
      DISN4S=SUMEXC(3,5,J)+DISNP*0.5
      DISN2D=SUMEXC(3,6,J)+DISNP*0.5
      PDISOP=SUMION(2,4,J)+SUMION(2,5,J)
      UVDISN=OTHPR1(1,J)
      PRHEP=OTHPR1(2,J)
      PO1DSR=OTHPR1(3,J)
      PO1SEL=SUMEXC(1,2,J)
      PN2PEL=SUMEXC(3,7,J)
      PO1DEL=SUMEXC(1,1,J)
C---- Calculate hv + NO ionization frequency from Lyman-a
      PLYNOP=OTHPR2(1,J)
C---- Calculate diss and ion rates for NOP - see also CODDN
C---- Calculate hv + NO dissociation frequency from hv + O2 dissoc rate
C---- Lyman-a NO+ freq used for NO dissoc by Lyman-a
      PSRNO=OTHPR2(2,J)
C
C----- Odd nitrogen zeroed to prevent convergence problems at low altitude
C----- or copied from the vertical grid
      IF(IVERT.GT.0) THEN
         N2D(J)=VN2D(J)
         N4S(J)=VN4S(J)
         NNO(J)=VNNO(J)
      ENDIF
      IF(Z(J).LT.100) NNO(J)=0.0
      IF(Z(J).LT.100) N4S(J)=0.0

      !-- Keep O+ density high to avoid convergence problem
      !..IF(Z(J).LT.900.AND.N(1,J).LT.1.0E-5) N(1,J)=N(1,J)+1.0E-5

C....... He+
       HEPLUS=N(4,J)
C
C....... XIONN is N+ calculated from the diffusion equation
       NPLUS(J)=XIONN(1,J)
C
C*****  obtain reaction rates from subr rats
      CALL RATS(J,TI(3,J),TI(1,J),TN(J),RTS)
C++++++++ multiply o+ + n2 rate by factor for vib n2 ++++
      IF(VC(J).LE.1) VC(J)=1.0
      RTS(3)=RTS(3)*VC(J)
      !.. RTS(99) is used in the RVN2PB file for N2+(v) calculation
      RTS(99)=RTS(42)*RTS(10)
C
C.......... Only recalculate minor ions if not already stored
      IF(I.EQ.1.OR.I.EQ.-1.OR.I.GT.2) GO TO 12
C
C........ use previously stored values  for [e] initially
      ZNE(J)=N(1,J)+N(2,J)+NOP(J)+O2P(J)+N2PLS(J)
      CALL CO2P(J,0,0,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),PO2P,
     >  DO2P,EUVION(2,7,J),PEPION(2,7,J),N(1,J),OP2D(J),N2PLS(J),
     >  NPLUS(J),N4S(J),NNO(J),OP2P(J))
C
C........... no+
      CALL CNOP(J,0,0,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),PNOP
     > ,DNOP,N(1,J),N2PLS(J),O2P(J),N4S(J),NNO(J),NPLUS(J),N2P(J)
     >  ,PLYNOP,VC(J),N2D(J),OP2D(J))
C
C.....calculation of [e] analytically for first guess in Newton
      B=N(1,J)+N(2,J)+N2PLS(J)
      C=PNOP/RTS(5)+PO2P/RTS(6)
      ZNE(J)=(B+DSQRT(B**2+4.0*C))/2.0
      ZNESAV=ZNE(J)
C
C....... call vibn2p to get effective n2+ + o->n2 + o+ rate
C...      IF(SZA(J).LT.1.57) CALL VIBN2P(J,ID,0,N2N(J),ON(J),ZNE(J),OP2P(J)
C...     > ,OP2D(J),RTS,....(3,1,J),...(3,2,J),....(3,3,J),O2N(J),NTOT,Z(J))
C
C.......... newton solution for electron density ZNE
C...... Loop through chemistry twice to evaluate dF/dh for Newton
      JITER=0
 9    DO 10 ITS=1,2
        IF(JITER.EQ.0) H=ZNE(J)*0.001
        IF(ITS.EQ.2) ZNE(J)=ZNE(J)+H
C............ o+(2p)
      CALL COP2P(0,J,0,0,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),OP2P(J)
     > ,TOP2PI,0.0D0,HEPLUS,N4S(J),NNO(J))
C............ o+(2d)
      CALL COP2D(J,0,0,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),OP2D(J)
     > ,EUVION(1,8,J),PEPION(1,8,J),OP2P(J),HEPLUS,N4S(J),NNO(J))
C
C........... n2+
      CALL CN2PLS(J,0,0,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,N2PLS(J),EUVION(3,1,J),EUVION(3,2,J),EUVION(3,3,J)
     > ,PEPION(3,1,J)+PAUION(3,1,J),PEPION(3,2,J)+PAUION(3,2,J)
     > ,PEPION(3,3,J)+PAUION(3,3,J),OP2D(J),OP2P(J),HEPLUS
     > ,NPLUS(J),NNO(J),N4S(J))
C
C........... o2+
      CALL CO2P(J,0,0,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),PO2P
     > ,O2P(J),EUVION(2,7,J),PEPION(2,7,J),N(1,J),OP2D(J),N2PLS(J),
     > NPLUS(J),N4S(J),NNO(J),OP2P(J))
C
C........... no+
      CALL CNOP(J,0,0,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),PNOP
     > ,NOP(J),N(1,J),N2PLS(J),O2P(J),N4S(J),NNO(J),NPLUS(J),N2P(J)
     > ,PLYNOP,VC(J),N2D(J),OP2D(J))
C
      FEX(ITS)=ZNE(J)-N(1,J)-N(2,J)-NOP(J)-O2P(J)-N2PLS(J)-
     >  OP2D(J)-OP2P(J)-NPLUS(J)-HEPLUS
 10   CONTINUE
C
      JITER=JITER+1
        IF(JITER.GT.22) THEN
         WRITE(6,'(A,I3)') '  ID=',ID
         WRITE(6,*) '  J   ALT      [e]    [O+]     [NO+]   [O2+]'
         WRITE(6,'(I3,F9.1,1P22E9.2)') J,Z(J),ZNE(J),N(1,J),
     >     NOP(J),O2P(J)
         CALL RUN_ERROR    !.. print ERROR warning in output files
         STOP
        ENDIF
      DEX=(FEX(2)-FEX(1))/H
      ZNE(J)=ZNE(J)-H-FEX(1)/DEX
      IF(ABS(FEX(1)/DEX/ZNE(J)).GT.1.0E-2) GO TO 9

 12   CONTINUE

C........n2(a3sigma)
      CALL CN2A(J,0,0,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,N2A,SUMEXC(3,1,J),SUMEXC(3,2,J),SUMEXC(3,3,J),SUMEXC(3,4,J))
C...... densities for O2 atmospheric bands
        CALL CO2(J,0,0,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,BOD3,O1D(J),O2SS,TN(J),O2B1,O2DEL,O2ADEL,O2ASIG)
C.......O(1S)
      CALL CO1S(J,0,0,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,O1S,O2P(J),N4S(J),PO1SEL,NPLUS(J),N2A,PO1DSR,O2ADEL,O2DEL)
C
C
C......... local equilibrium densities for metastables and odd nitrogen 
C
C....... N(2P) density
      CALL CN2P(J,0,0,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),PROD(5)
     > ,LOSS(5),N2P(J),PN2PEL,UVDISN,O2P(J),NNO(J),N2PLS(J))
      IF(LOSS(5).GT.0) N2P(J)=PROD(5)/LOSS(5)
C
C...... N(2D) density
      CALL CN2D(J,0,0,Z(J),RTS,ON(J),O2N(J),N2N(J),NOP(J),ZNE(J),PROD(1)
     > ,LOSS(1),N2PLS(J),DISN2D,UVDISN,NPLUS(J),N2P(J)
     > ,N2D(J),N(1,J),NNO(J),N2A)
      N2D(J)=PROD(1)/LOSS(1)
      IF(IVERT.GT.0) N2D(J)=VN2D(J)
C
C......... NO density
      CALL CNO(J,0,0,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),PROD(3),LOSS(3)
     >  ,N2D(J),N4S(J),N2P(J),NNO(J),O2P(J),N(1,J),PSRNO,PLYNOP,N2A,
     >  NPLUS(J))
      IF(LOSS(3).GT.0) NNO(J)=PROD(3)/LOSS(3)
      IF(IVERT.GT.0) NNO(J)=VNNO(J)
C
C--- N4S is taken from MSIS in CMINOR if not calculated on vertical grid
      IF(IVERT.GT.0) N4S(J)=VN4S(J)
C
C........ O(1D) density
      CALL CO1D(0,J,0,0,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J),PROD(4)
     > ,LOSS(4),O2P(J),N2D(J),PO1DEL,PO1DSR,O1S,NPLUS(J),O1D(J))
      IF(LOSS(4).GT.0) O1D(J)=PROD(4)/LOSS(4)
C
      FD(8)=N2D(J)*RTS(8)*ZNE(J)
      FD(1)=NOP(J)
      FD(2)=O2P(J)
      FD(3)=N2PLS(J)
      FD(4)=OP2P(J)
      FD(5)=OP2D(J)
      FD(6)=OP2D(J)*RTS(12)*ZNE(J)
      FD(7)=NPLUS(J)
C........ calc o+ prod. from minor ions
      P1=OP2P(J)*(RTS(26)*ON(J)+RTS(14)*ZNE(J)+0.047)
      P2=OP2D(J)*(RTS(12)*ZNE(J)+RTS(28)*ON(J))
      P3=N2PLS(J)*RTS(99)*ON(J)
      P4=NPLUS(J)*ON(J)*RTS(31)
      FD(9)=P1+P2+P3+P4
      IF(I.LT.2) RETURN
C
C
C...... this next section for printing densities and production rates
      DISN2S=SUMION(3,6,J)
      P3X12=SUMEXC(3,12,J)
      P3X11=SUMEXC(3,11,J)
      P3X2=SUMEXC(3,2,J)
      P3X3=SUMEXC(3,3,J)
      P1X3=SUMEXC(1,3,J)
      P1X4=SUMEXC(1,4,J)
      P1X5=SUMEXC(1,5,J)
      P1X6=SUMEXC(1,6,J)
      P1X7=SUMEXC(1,7,J)
      P1X8=SUMEXC(1,8,J)
      P1X9=SUMEXC(1,9,J)
C
      IF(I.EQ.2) CALL CION(J,ID,JP,Z(J),N(1,J),O2P(J),NOP(J),N2PLS(J)
     > ,NPLUS(J),OP2D(J),OP2P(J),ZNE(J),TI(3,J),TI(1,J),VC(J),HEPLUS
     > ,N(2,J),NNO(J),TN(J))
C
      IF(I.EQ.3) CALL CNEUT(J,ID,JP,Z(J),O1S,N2A,N2D(J),O1D(J),NNO(J),
     > N4S(J),N2P(J))
C
      IF(I.EQ.4) THEN
        CALL CEMM(J,ID,JP,Z(J),RTS,O1D(J),N2D(J),OP2D(J),OP2P(J),O1S
     >   ,DISNP,ZNE(J),ON(J),N2N(J),N4S(J),P3X12,P3X11,N2A
     >   ,P3X2,P3X3,DISN2S)
      ENDIF
C
      IF(I.EQ.5) CALL CO2EM(J,ID,JP,Z(J),O2B1,O2DEL,O2SS,O2ADEL,O2ASIG)
C
      IF(I.EQ.6) CALL CN2D(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),NOP(J)
     > ,ZNE(J),PROD(1),LOSS(1),N2PLS(J),DISN2D,UVDISN,NPLUS(J),N2P(J)
     > ,N2D(J),N(1,J),NNO(J),N2A)
C
C------- O(1D) production and loss rates
      IF(I.EQ.7) CALL CO1D(0,J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J)
     > ,ZNE(J),PROD(4),LOSS(4),O2P(J),N2D(J),PO1DEL,PO1DSR,O1S
     > ,NPLUS(J),O1D(J))
C
C------- O(1D) volume emission rates
      IF(I.EQ.36) CALL CO1D(1,J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J)
     > ,ZNE(J),PROD(4),LOSS(4),O2P(J),N2D(J),PO1DEL,PO1DSR,O1S
     > ,NPLUS(J),O1D(J))
C
      IF(I.EQ.8) CALL CNO(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,PROD(3),LOSS(3),N2D(J),N4S(J),N2P(J),NNO(J),O2P(J),N(1,J)
     >  ,PSRNO,PLYNOP,N2A,NPLUS(J))
C
      IF(I.EQ.9) CALL CN4S(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,PROD(2),LOSS(2),N4S(J),DISN4S,N2D(J),N2P(J),N(1,J),N2PLS(J)
     > ,UVDISN,NOP(J),NPLUS(J),NNO(J),O2P(J),PSRNO,VC(J))
C
      IF(I.EQ.10) CALL CO1S(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,O1S,O2P(J),N4S(J),PO1SEL,NPLUS(J),N2A,PO1DSR,O2ADEL,O2DEL)
C
      IF(I.EQ.11) CALL CN2PLS(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J)
     > ,ZNE(J),N2PLS(J),EUVION(3,1,J),EUVION(3,2,J),EUVION(3,3,J)
     > ,PEPION(3,1,J)+PAUION(3,1,J),PEPION(3,2,J)+PAUION(3,2,J)
     > ,PEPION(3,3,J)+PAUION(3,3,J),OP2D(J),OP2P(J),HEPLUS
     > ,NPLUS(J),NNO(J),N4S(J))

      IF(I.EQ.12) CALL CNOP(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,PNOP,NOP(J),N(1,J),N2PLS(J),O2P(J),N4S(J),NNO(J),NPLUS(J)
     > ,N2P(J),PLYNOP,VC(J),N2D(J),OP2D(J))
C
      IF(I.EQ.13) CALL CO2P(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,PO2P,O2P(J),EUVION(2,7,J),PEPION(2,7,J),N(1,J),OP2D(J),N2PLS(J),
     > NPLUS(J),N4S(J),NNO(J),OP2P(J))
C
      IF(I.EQ.14) THEN
        EUVP4S=EUVION(1,7,J)
        PEOP4S=PEPION(1,7,J)+PAUION(1,7,J)
        CALL COP4S(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,N(1,J),EUVP4S,OP2D(J),OP2P(J),PEOP4S,PDISOP,N2PLS(J)
     > ,N2D(J),NNO(J),VC(J),HEPLUS,NPLUS(J))
      ENDIF
C
      IF(I.EQ.15) CALL COP2D(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     >,OP2D(J),EUVION(1,8,J),PEPION(1,8,J),OP2P(J),HEPLUS,N4S(J),NNO(J))
C
      IF(I.EQ.16) THEN
        PEOP2P=PEPION(1,9,J)+PAUION(1,9,J)
        CALL COP2P(1,J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     >   ,OP2P(J),TOP2PI,PEOP2P,HEPLUS,N4S(J),NNO(J))
      ENDIF
      !.. Get 7320 volume emission rates
      IF(I.EQ.22) THEN
        PEOP2P=PEPION(1,9,J)+PAUION(1,9,J)
        CALL COP2P(1,J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     >   ,OP2P(J),TOP2PI,PEOP2P,HEPLUS,N4S(J),NNO(J))
      ENDIF
C
      IF(I.EQ.17) CALL CNPLS(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,DISNP,PEDISNP,NPLUS(J),N(1,J),N2D(J),OP2P(J),HEPLUS,OTHPR2(3,J)
     > ,O2P(J),N4S(J),OP2D(J),N2PLS(J),NNO(J))
C
      IF(I.EQ.18) CALL CN2A(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,N2A,SUMEXC(3,1,J),SUMEXC(3,2,J),SUMEXC(3,3,J),SUMEXC(3,4,J))
C
      IF(I.EQ.19) CALL CN2P(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,PROD(5),LOSS(5),N2P(J),PN2PEL,UVDISN,O2P(J),NNO(J),N2PLS(J))
C
      IF(I.EQ.20) CALL CHEP(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,HEPLUS,PRHEP,HE,NNO(J),N(1,J),N(2,J),HN(J))
C
      IF(I.EQ.21) CALL GRAV_WAV(J,ID,JP,Z(J),ON(J),O2N(J),N2N(J),HE
     > ,HN(J),TN(J))
C
      !...IF(I.EQ.22) CALL CRT2(J,ID,JP,Z(J),RTS) !.. reaction rates unused
C
      IF(I.EQ.23) CALL CMSIS(J,ID,JP,Z(J),ON(J),O2N(J),N2N(J),HE
     > ,HN(J),TN(J),O1S,N2A,N2D(J),O1D(J),NNO(J),N4S(J),N2P(J))
C
      IF(I.EQ.24) CALL C1200(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,P3X12,P3X11,N4S(J),DISNP)
C
      IF(I.EQ.25) CALL CO2(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     > ,BOD3,O1D(J),O2SS,TN(J),O2B1,O2DEL,O2ADEL,O2ASIG)
C
      IF(I.EQ.26) THEN
        PN2PX=SUMION(3,1,J)
        PN2PA=SUMION(3,2,J)
        PN2PB=SUMION(3,3,J)
        CALL CN2PV(J,ID,JP,Z(J),RTS,ON(J),O2N(J),N2N(J),ZNE(J)
     >   ,OP2P(J),OP2D(J),PN2PX,PN2PA,PN2PB)
      ENDIF
C
      IF(I.EQ.27) CALL PRDPRB(J,ID,JP,Z(J))
C
      !.. Electron impact excitation rates for atomic oxygen.Modified 2020-08-02 
      !.. to include O1D to calculate 6300 emission rates
      IF(I.EQ.28) CALL CEMOX(J,ID,JP,Z(J),RTS,PO1DEL,PO1SEL,P1X3,P1X4,
     > P1X5,P1X6,P1X7,P1X8,P1X9,OTHPR2(4,J),ZNE(J),N(1,J),O1D(J),O1S)
C
C...... This part to print O2+ production rates DEACTIVATED (use I=28)
      IF(I.EQ.-28) THEN
       IF(JP.EQ.1) WRITE(ID,*) ' ' 
       IF(JP.EQ.1) WRITE(ID,*) '  ALT    EUV O2+'
        WRITE(ID,777) Z(J),(EUVION(2,K,J),PEPION(2,K,J),K=1,5)
      ENDIF
 777    FORMAT(F6.1,1P,22E8.1)
      RETURN
      END
